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U«IGOR 


ABSTRACT 


In  the  first  phase  of  the  project  to  develop  an  articulated  human  biomechanical  modeling 
toolbox  (AHBM),  efforts  were  made  to  review  the  current  status  of  modeling  techniques  and  to 
develop  rigid  body  formulations  suitable  for  human  biomechanical  modeling.  The  developed 
version-one  of  the  toolbox,  AHBM  VI,  consists  of  kinematics,  inverse  and  forward  dynamics 
algorithms,  data  conversion  routines,  graphical  algorithms,  and  many  utility  routines  for 
mathematical  calculation  and  file  operation.  Both  inverse  and  forward  sample  models,  such  as  a  3D 
lower  extremity  model,  a  3D  whole  body  human  model,  and  a  human  head-neck  model,  were 
developed.  AHBM  VI  was  also  used  for  two  application  projects:  upgrading  gait  analysis  software 
for  USARIEM;  and  developing  a  new  method  to  understand  the  airbag  external  load  behavior  and 
bag-occupant  interaction  for  Department  of  Transportation. 

The  following  specific  tasks  were  completed  in  the  first  phase 

♦  Review  of  current  techniques  for  human  biomechanical  modeling 

♦  Developing  three-dimensional  kinematics  algorithms 

♦  Developing  forward  dynamics  formulations  for  multibody  constrained  systems 

♦  Developing  commonly  used  algorithms  related  to  inverse  dynamics  analysis 

♦  Developing  data  structure  and  data  conversion  routines 

♦  Developing  some  graphical  algorithms  for  visualizing  data  and  model 

♦  Developing  utility  routines  for  mathematical  calculation  and  file  operation 

♦  Developing  example  problems  for  both  inverse  and  forward  dynamics  analysis 

In  the  phase  two  of  the  project,  efforts  will  be  focused  on  reviewing  and  developing  muscle 
models.  The  muscle  models  wiU  be  included  in  AHBM  V2.  This  will  expand  the  toolbox’s 
capabilities  to  handle  human  biomechanical  modeling  where  muscle  activities  have  to  be  accounted 
for.  Specific  applications  may  include  solving  muscle  load  sharing  in  inverse  dynamics  analysis  and 
understanding  overuse  injuries  related  to  muscle  activities. 

This  report  is  separated  into  two  parts.  The  first  part  consists  of  an  overview  of  the  toolbox, 
rigid  body  formulations,  and  example  models  and  applications.  The  second  part  provides  detailed 
description  of  individual  routines  in  the  toolbox. 
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CHAPTER  1 
OVERVIEW 


1.1  INTRODUCTION 

This  report  documents  the  progresses  on  developing  a  human  biomechanical  modeling 
toolbox  during  the  first  phase  of  the  research  project.  It  is  separated  into  two  parts.  This  first  part 
consists  of  an  overview  of  the  toolbox,  rigid  body  formulations,  and  example  models  and 
applications.  A  separate  second  part  provides  detailed  description  of  individual  routines  in  the 
toolbox. 

This  first  part  is  divided  into  five  chapters.  The  first  chapter  gives  an  overview  of  the  current 
state  of  the  analysis,  modeling  and  simulation  of  human  biomechanical  system.  The  strategy  and 
plan  for  developing  a  toolbox  is  also  given.  Chapter  two  gives  the  formulations  for  kinematics 
calculation  with  an  emphasis  on  three-dimensional  rotational  kinematics.  Chapter  three  develops 
forward  dynamics  formulations  for  multibody  systems  that  are  suitable  for  forward  dynamic  simula¬ 
tion  of  human  systems.  Chapter  four  describes  the  topics  commonly  associated  to  the  inverse 
dynamic  analysis  of  human  biomechanical  systems.  Chapter  five  provides  some  example  models 
and  applications  of  the  current  version  of  the  toolbox. 

1.2  BIOMECHANICAL  ANALYSIS,  MODELING  AND  SIMULATION 
1.2.1  Need  for  Biomechanical  Modeling 

The  analysis,  modeling  and  simulation  of  human  system  have  become  a  more  and  more 
important  research  area.  This  is  driven  by  our  intrinsic  curiosity  of  understanding  the  fundamental 
mechanisms  as  well  as  by  the  needs  from  both  the  civil  and  the  mhitaiy.  For  the  military,  it  always 
faces  the  challenge  of  evaluating  and  improving  soldiers’  performance  under  various  environmental 
and  operational  conditions,  as  well  as  reducing  associated  injuries.  Biomechanical  researchers  in  the 
military  have  to  deal  with  issues  such  as  developing  better  equipment  to  enhance  soldiers’ 
performance,  design  better  training  schedule  to  increase  soldiers’  strength  without  incurring  injury, 
designing  better  ways  of  monitoring,  assessing  and  improving  the  medical  status  of  warfighters 
during  training  and  combats.  A  very  important  part  of  these  research  issues  is  the  analysis  of 
mechanical  loads  and  energetic  requirements  of  the  human  body  during  various  movements. 

Traditionally,  most  biomechanical  research  involves  only  experimentation  where  the 
statistics  of  laboratory  measurements  and/or  field  observations  is  the  primary  research  result.  The 
experimental  approach,  however,  suffers  several  limitations.  First,  human  body  is  a  complex,  highly 
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nonlinear  system.  Significant  variations  of  the  measurements  of  human  responses  are  expected, 
especially  among  different  individuals.  And  in  many  cases  the  data  cannot  be  repeatable  with  good 
accuracy.  Therefore,  the  applicability  of  many  experiment  results  has  to  be  examined  with  caution. 
Second,  even  when  statistically  significant  relationships  can  be  established  from  experiments,  the 
underlying  mechanisms  are  usually  too  complex  to  be  revealed.  Therefore,  most  experiment  results 
can  only  be  used  to  answer  specific  questions.  They  usually  do  not  provide  insight  into  the  problems 
themselves.  Third,  experimentally  work  are  often  expensive  and  difficult  to  control,  taking  into 
account  that  a  large  number  of  people  over  a  long  period  of  time  are  usually  required  to  get 
statistically  significant  results.  Finally  in  many  cases  it  is  ethically  impossible  to  perform  experiments 
on  human  beings,  especially  when  injuries  may  be  involved  in  the  tests. 

WhUe  experimental  work  stiU  remains  and  will  always  be  an  important  part  of 
biomechanical  research,  more  and  more  efforts  are  being  put  into  the  use  of  modeling  and 
simulation  in  solving  biomechanical  problems  associated  with  human  systems.  In  modeling,  the 
human  system  is  represented  by  sets  of  mathematical  relationships  including  some  modeling 
parameters.  By  varying  model  parameters,  numerical  simulation  can  be  performed.  Compared  with 
experimental  work,  numerical  simulations  are  usually  faster,  cheaper  and  easier  to  control.  The 
models  can  be  used  for  various  situations  and  the  results  are  more  insightful. 

Whether  a  model  can  be  successfully  used  for  a  practical  biomechanical  problem  depends  on 
how  accurate  the  model  imitates  the  real  problem,  whether  there  are  enough  mathematical  tools  and 
computational  power  to  implement  the  model,  and  whether  there  are  enough  experimental  data  to 
verify  the  model  results.  A  combination  of  experimentation,  biomechanical  analysis  and 
biomechanical  modeling  and  simulation  is  usually  the  best  approach  to  tackle  complex  real  world 
problems.  The  following  sections  will  brief  summarize  the  current  status  of  biomechanical  analysis, 
modeling  and  simulation  of  human  movement. 

1.2.2  Model  Sub-systems 

Human  as  a  living  biomechanical  system  is  a  complex  and  integrated  neural-muscular- 
skeletal  system  combing  motion  generation  components,  adaptive  control,  reflexes,  self-analysis  and 
learning  (Barnes,  Oggero  et  al.  1997).  Any  voluntary  human  movements  are  complex  motions 
characterized  by  the  presence  of  the  so-called  controlling  and  compensating  parts  as  the  result  of  the 
coordination  of  neural  control,  muscle  activation,  and  skeletal  motion  generation.  The  central  neural 
system  sets  the  goal  to  be  achieved.  The  neuromuscular  system  correspondingly  control  the  pattern 
of  the  muscle  activation  generate  tension.  The  muscle  forces  drive  the  skeletal  system  to  generate  or 
adjust  the  motion. 

In  modeling,  the  different  levels  of  the  human  biomechanical  system  are  represented  by 
different  levels  of  mathematical  systems. 
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Mechanical  Subsystem 

Rigid  body  representation  of  the  human  skeletal  system  is  commonly  used  where  rigid 
bodies  with  masses  are  connected  at  articulations.  The  governing  equation  of  motions  can  be 
obtained  from  the  Newton-Euler  formations,  Lagrange  formulations  or  Kane’s  method.  The  resulted 
sets  of  equations  include  kinematic  data  (velocities  and  accelerations),  kinetic  data  (forces  and 
moments)  and  inertial  properties  of  body  segments  (mass  and  moment  of  inertia). 

Deformable  bodies  are  have  also  been  used  in  some  models.  For  example,  long  bones  can  be 
represented  as  deformable  beams.  Finite  element  models  and  other  models  based  on  continuum 
mechanics  are  also  used  to  model  the  mechanical  subsystem.  Simulation  of  specific  structures,  such 
as  the  foot,  benefits  from  the  use  of  deformable  body  representation,  but  due  to  the  limitations  of 
modeling  complexity,  rigid  body  formulation  are  mostly  used  in  biomechanical  analysis. 

Actuator  Subsystem 

Muscles,  which  are  the  only  actuators  that  can  generate  forces,  are  modeled  by  actuator 
subsystem  governed  by  muscle  dynamics  and  muscle  metabolics.  The  formulations  for  the  actuator 
subsystem  include  models  of  various  complexities.  Some  involves  only  the  passive  elements  and  no 
active  actuators  are  present.  More  sophisticated  muscle  models  involve  both  the  passive  and  the 
active  elements.  They  also  include  the  muscle  insertion  points  in  order  to  determine  the  lines  of 
action  of  the  muscle  forces.  These  components  obviously  become  highly  integrated  with  the 
geometry  selected  for  the  physical  representation  of  the  body  segments.  The  difficulties  in 
implementing  these  models  include  the  modeling  of  multiarticular  muscles  that  span  over  more  than 
one  joint  and  the  redundancies  presented  by  the  multiple  muscles  on  single  joints.  Even  more 
complicated  models  of  individual  muscles  have  also  been  postulated  to  understand  muscle  function. 
Individual  muscle  spindles  and  even  individual  motor  units  and  their  recruitment  and 
neuromuscular  phenomena  have  been  modeled.  However,  currently,  these  models  only  serve  to 
provide  insight  into  how  exactly  muscles  function. 

Control  Subsystem 

The  control  subsystem  with  a  variety  of  control  schemes  simulates  the  function  of  the  neural 
system.  Control  strategies  of  various  level  of  complexity  have  been  developed.  The  simplest  strategy 
is  to  have  no  actuators  to  control.  Among  controlled  models,  open-looped  control  and  closed-loop 
control  has  been  studied  extensively.  More  adaptive  approaches  using  neural  network  and  feed 
forward  control  are  also  under  intensive  research. 

1.2.3  Types  of  Biomechanical  Analysis  and  Modeling 

Depends  on  the  goal  of  a  specific  biomechanical  problem,  various  types  of  analysis  and 
modeling  techniques  can  be  adopted. 
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Inverse  Dynamics  Approach 

Inverse  d)rnamics  approach  involves  the  mechanical  subsystem  and  possibly  the  actuator 
subsystem.  At  the  mechanical  level,  it  calculates  the  joint  forces  and  moments  from  segmental 
inertial  properties,  measured  kinematics  data  and  external  forces. 

The  necessary  kinematic  data  can  be  measured  from  accelerometers  mounted  on  body 
segments.  The  acceleration  data  is  then  integrated  to  get  the  required  velocities  and  displacements. 
More  frequently,  displacement  data  are  obtained  from  a  motion  capture  system  (video-based  or 
optoelectronic).  The  displacement  data  are  then  differentiated  to  obtain  the  velocities  and 
accelerations. 

Inertia  properties  of  body  segments  can  be  obtained  through  direct  measurement  on  the 
subject.  There  are  also  regression  equations  based  on  a  number  of  anthropometric  measurements 
such  as  body  v^^eight,  stature  and  specific  geometric  dimensions  of  individual  segments. 

The  external  forces  in  the  biomechanical  analysis  of  human  systems  are  usually  ground 
reaction  forces.  Force  plates  can  measure  the  ground  reaction  forces.  The  measured  ground  reaction 
forces  usually  have  to  be  processed  (filtered,  synchronized  with  kinematics  data)  before  can  be  used 
for  inverse  dynamics  analysis. 

At  the  mechanical  level,  the  inverse  dynamics  analysis  provides  the  Joint  forces  and 
moments.  Some  energetics  quantities  such  as  internal  work,  extemed  work,  Joint  power  flow  may 
also  be  obtained.  These  kinetics  and  energetics  quantities  are  used  in  many  researches  to  serve  the 
purposes  such  as  analyzing  both  the  healthy  and  pathological  gait,  assisting  in  the  diagnosis  of  the 
underlying  pathologies  of  abnormal  gait  pattern,  helping  the  evaluation  and  design  of  soldier 
equipment,  etc.  In  many  of  these  researches,  statistics  analysis  of  test  results;  characterization  and 
pattern  analysis  of  the  time  history  of  the  results  have  to  be  used. 

After  the  Joint  forces  and  moments  are  obtained,  they  can  be  used  to  estimate  muscle  forces 
and  Joint  reaction  forces  from  passive  structures  (ligaments,  tendon).  Since  a  large  number  of 
muscles  spanning  each  Joint  in  the  human  body,  the  estimation  of  muscle  forces  is  an  indeterminate 
problem  (Collins  1995).  Techniques  to  solve  this  problem  are  based  on  either  grouping  muscles  with 
similar  function  thus  eliminating  redundancy  or  applying  optimization  criteria  to  solve  the  muscle 
force  distribution  (Crowninshield  1978;  Hardt  1978;  Patriarco,  Mann  et  al.  1981;  Vaughan,  Hay  et 
al.  1982;  Cholewicki  and  McGill  1994).  The  first  approach  leads  to  oversimplification  and 
unsatisfactory  results.  The  second  approach  is  usually  called  inverse  optimization  or  static 
optimization.  The  inverse  optimization  are  numerical  efficient  and  has  been  successfully  applied  but 
it  also  suffers  from  certain  drawbacks.  First  the  optimization  criteria  has  not  been  completely 
understood.  Second,  the  inverse  optimization  cannot  guarantee  the  continuity  of  the  solution.  Third, 
muscle  dynamics  is  not  reflected  in  inverse  optimization.  Recently,  some  work  (Happee  1994)  has 
been  done  to  include  muscle  dynamics  as  constraints  in  inverse  optimization,  which  has  been  shovm 
to  work  well  to  solve  complex  shoulder  mechanism  (Happee  and  Van  der  Helm  1995).  Finally  the 


lack  of  reliable  validation  procedures  cilso  limits  the  application  of  inverse  optimization.  Currently, 
electromyogram  signals,  which  describe  the  input  into  the  muscular  system  is  usually  recorded  to  see 
if  it  matches  with  the  muscle  force  pattern  calculated  from  inverse  dynamics  optimization. 

Forward  Dynamics  Approach 

While  inverse  d5mamics  approach  tries  to  calculate  joint  forces  and  moments  and  thus 
muscle  forces  from  measured  kinematics  data,  forward  dynamics  approach  attempts  to  predict  the 
motion  from  neural  or  mechanical  inputs.  A  forward  simulation  may  involve  aU  the  mechanical, 
actuator  and  the  control  systems. 

When  no  actuator  is  present,  only  the  initial  conditions  are  required  and  the  simulation 
predicts  the  motion  trajectories.  This  involves  only  the  mechanical  subsystem.  The  approach  has 
been  widely  used  in  simulating  some  traumatic  events,  such  as  car  crash  where  the  event  happen  in  a 
very  short  period  time  and  the  effects  of  neural  control  and  muscle  activities  can  be  neglected. 

As  a  control  model,  the  easiest  one  to  implement  is  an  open-loop  controller  with  a  torque 
that  is  only  a  function  of  time.  This  removes  all  the  complexities  associated  with  computing 
kinematics  of  muscle  insertion  points.  How  the  joint  torque  changes  with  respect  to  time  depends 
mainly  on  the  goal  of  the  simulation.  One  way  is  to  solve  the  inverse  dynamic  problem  using  data 
collected  during  an  experiment  in  which  a  subject  performs  the  same  actions  that  the  simulation  is 
trying  to  reproduce.  Another  way  to  obtain  the  torque  is  to  impose  an  optimization  criterion.  The 
use  of  such  a  simplified  actuator  and  control  model  cannot  compensate  for  external  perturbations 
and  provides  no  means  for  actively  maintaining  balance  when  equilibrium  conditions  are  desired. 

When  force  and  torque  application  is  also  a  function  of  position,  velocity  or  other 
generalized  variables  in  the  model,  the  control  become  closed  loop.  Muscle  and  spinal  level  reflex  is 
an  example  of  closed-loop  control.  Simulation  of  higher  level  controls,  such  as  cerebellum  and  brain, 
becomes  more  challenging.  More  adaptive  approaches,  such  as  neural  networks  and  artificial 
intelligence  techniques,  can  be  applied  to  model  these  interactions  and  relate  other  biomechanically 
relevant  parameters  to  control. 

Feedforward  concepts  are  also  adopted  in  the  control  theory.  A  feedforward  controller 
(brain)  predicts  the  future  motion  of  the  controlled  system  (the  body)  and  compares  it  to  the  current 
body  position.  When  disaepancy  is  noted,  appropriate  corrective  signals  are  issued.  A  key  point  of 
feedforward  control  is  that  it  is  a  dynamic  process  and  this  can  be  adapted  to  explain  some  processes 
of  motor  learning  and  performance  enhancement.  As  increasing  possible  perturbations  have  been 
experienced  during  the  repetitions  or  training,  patterned  responses  are  learned.  The  feedforward 
control  requires  a  very  accurate  and  rapid  motion  analysis.  Currently,  most  simulations  of  human 
system  do  not  incorporate  movement  analysis  capability.  However,  further  development  of 
feedforward  concept  will  lead  to  more  useful  simulations. 


5 


Energetics 

The  energetic  aspect  human  biomechanical  system  deals  with  the  mechanical  power  as  well 
as  the  metabolic  cost  associated  with  human  movements.  The  metabolic  cost  can  be  separated  into 
different  terms  according  to  their  sources.  The  activation  heat  accounts  for  the  energy  used  to 
activate  muscles.  The  maintenance  heat  is  used  to  maintenance  the  muscle  forces.  The  shorting  heat 
is  related  to  the  extra  heat  produced  as  a  consequence  of  the  shortening  of  muscle.  The  mechanical 
work  is  the  product  of  muscle  force  doing  work.  And  the  dissipation  of  energy  in  passive  structures 
also  contributes  to  the  total  metabolic  costs.  Although  research  work  has  led  to  some  empirical 
relations,  the  actual  metabolic  cost  is  usually  determined  from  the  measured  oxygen  consumption, 
which  is  related  to  the  metabolic  cost. 

Different  methods  are  available  to  estimate  the  mechanical  power  during  human  movement. 
Some  methods  calculate  the  power  based  on  external  work  necessary  to  move  the  center  of  mass  of  a 
human  body.  However,  the  actual  dynamics  of  each  segment  is  lost  in  these  methods.  Some  other 
methods  provide  the  mechanical  power  due  to  the  movement  of  each  segment.  These  models  cannot 
account  for  the  synergy  of  muscles  over  a  single  joint.  Although  it  is  recognized  that  the  actual 
mechanical  work  and  thus  the  metabolic  cost  has  to  come  from  the  muscles,  currently  it  has  not 
been  very  successful  in  calculating  muscle  power  due  to  the  difficulties  in  determining  muscle  forces. 

1 .3  DEVELOPING  A  BIOMECHANICAL  MODELING  TOOLBOX 

Human  biomechanical  systems  are  highly  non-linear  systems  with  a  high  number  of 
interconnected  and  interacting  elements.  There  are  multilevel  and  multi-loop  controllers  involved 
and  the  systems  are  characterized  by  highly  functional  adaptations.  Summary  of  the  available 
theoretical  and  experimental  investigations  on  the  human  biomechanical  systems  convinces  one  that 
a  complex  system  approach  has  to  be  adopted.  System  theory  implies  integral  analysis  of  the 
working  system  instead  of  studying  of  separate  phenomenon.  This  requires  the  biomechanical 
models  be  simple  to  use,  goal  oriented,  reliable,  complete  and  adaptive.  The  basis  of  this  is  a 
complex  of  mathematical  tools  that  can  be  assembled  to  describe  the  processes  in  real  human 
biomechanical  systems  under  investigation. 

On  the  other  hand,  although  biomechanical  problems  are  versatile  in  nature.  Some  basic 
mathematical  techniques  are  shared.  Providing  a  set  of  these  mathematical  tools  will  facilitate  and 
standardize  the  solution  of  these  complex  problems  without  reworking  on  the  mathematical  details 
which  is  usually  very  time  consuming  and  error  prone. 

Although  a  large  amount  of  commercial  and  research  software  have  been  developed  or  used 
in  biomechanical  analysis  of  human  systems.  Most  of  them  are  targeted  at  solving  specific  classes  of 
problems.  And  none  of  them  provides  such  a  basic  set  of  mathematical  tools.  It  is  the  object  of  this 
research  project  to  review,  collect,  organize  and  develop  these  fundamental  mathematical  tools 
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serving  human  biomechanical  analysis,  modeling  and  simulation.  A  human  modeling  toolbox  will 
be  developed  with  the  following  features 

1 .  Flexible  and  powerful:  users  wiU  be  provided  with  the  fundamental  and  powerful  tools  that 
can  be  used  to  solve  very  general  problems.  The  toolbox  will  be  organized  in  such  a  way  that 
models  can  easily  be  assembled  for  problems  under  investigation.  The  results  of  simulation 
can  be  visualized.  Attention  will  also  be  paid  on  the  development  of  standard,  convenient 
information  input  and  editing  facilities 

2.  Reliable:  the  mathematical  and  numerical  algorithms  will  be  rigorously  tested  and  well 
studied,  thus  avoiding  unnecessary  errors  in  mathematical  modeling  which  is  the  major 
source  of  error  in  biomechanical  modeling 

3.  Efficient:  complicated  models  can  be  assembled  from  standardized  functions.  Therefore, 
users  are  spared  from  basic  programming  and  can  concentrate  on  the  modeling  itself.  In 
addition,  higher-level  routines  specifically  developed  for  common  problems  will  also  be 
provided 

4.  Consistent:  developing  model  in  a  systematic  way  wiU  make  the  comparison  of  results  easier. 
Also,  benchmark  data  and  problems  wUl  also  be  included  in  the  toolbox  for  easy  calibration 
of  model  results 

5.  Open:  the  toolbox  wUl  allow  the  addition  of  new  algorithms  and  models.  This  is  crucial  for 
the  fast-growing  biomechanical  modeling  research  area.  The  open  structure  of  the  toolbox 
ensures  that  it  will  become  more  and  more  powerful  in  real  applications. 

6.  Easy  sharing  of  data:  data  wiU  be  stored  in  a  systematic  and  easUy  retrievable  way. 

In  order  to  satisfy  these  requirements,  the  foUowing  components  will  be  included  in  the  toolbox 

1.  Rigid  dynamics  routines:  including  kinematics,  inverse  dynamics,  and  forward  dynamics 
formulations  of  multibody  systems  with  closed  loop  and  constraints. 

2.  Muscle  modeling  routines:  including  various  models  for  the  passive  and  active  elements  of 
muscle  as  weU  as  a  collection  of  muscle  parameters. 

3.  Control  schemes  and  optimization  routines:  including  inverse  (static)  optimization, 
d)mamics  optimization,  open-looped  and  close-looped  control  using  both  joint  torque  and 
muscle  activation  as  control  parameters. 

4.  Data  analysis  routines:  including  the  smoothing,  fitting  and  filtering  of  time  history  data; 
statistical  analysis;  pattern  recognition  of  data. 

5.  Graphics  routines:  including  plotting  routines  for  data  viewing,  animation  routines  as  well  as 
routines  for  developing  graphic  user  interfaces 

6.  Other  utility  routines:  including  but  not  limited  to  routines  for  basic  matrix  manipulation, 
file  input  and  output,  data  format  conversion  etc. 
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7.  Database:  including  collected  model  parameters  and  simulation  results  as  well  as  routines  of 

managing  the  data 

Matlab  (Mathworks  2000)  is  an  integrated  technical  environment  that  combines  numeric 
computation,  advanced  graphics  and  visualization,  and  a  high-level  language.  This  makes  it  suitable 
for  the  analysis,  modeling  and  simulation  of  human  systems.  The  human  modeling  toolbox  will  be 
developed  in  Matlab.  However,  the  use  of  the  toolbox  will  not  be  limited  in  Matlab  environment. 
Programs  developed  in  Matiab  can  be  compiled  into  standalone  applications  as  well  as  web  based 
application. 

The  following  figure  shows  the  schematics  of  the  toolbox.  Key  to  the  toolbox  is  the  lower 
level  kernel  routines  that  perform  the  above  mentioned  functions.  Higher  lever  modeling  routines 
are  assembled  from  kernel  functions  for  commonly  encountered  modeling  problems.  Benchmark 
problems  will  be  developed  based  on  the  modeling  routines.  The  shaded  boxes  in  Figure  1-1, 
including  user  models  and  applications,  are  not  the  parts  of  the  toolbox. 

The  toolbox  will  be  developed  in  three  phases.  The  first  phase  will  focus  on  the  rigid  body 
formulation.  The  second  and  the  third  phases  deal  with  muscle  modeling  and  control  schemes  of 
different  levels  of  complexity.  The  development  of  graphics  and  additional  utility  routines  will 
spread  over  aU  the  three  phases. 
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Application  _ Algorithm  Deveiopment 


KERNEL  ROUTINES 


Rigid  Dynamics  Solver 

♦  Kinematics  routines 

♦  Inverse  dynamics  routines 

♦  Forward  dynamics  routines 


Muscle  dynamics 

♦  Passive  elements 

♦  Active  elements 

♦  Muscle  library 


Control  and  Optimization 

♦  Static  optimization 

♦  Dynamics  optimization 

♦  Open/close-looped  control 
using  joint  torques 

♦  Open/close-looped  control 
with  muscle  activation 

♦  Adaptive  control,  feed¬ 
forward  control,  etc. 


Data  Analysis 

♦  Smoothing,  fitting  filtering 

♦  Statistics  analysis 

♦  Pattern  recognition 


Graphics 

♦  2D,  3D  graphics 

♦  Animation 

♦  User  interface  tools 


DataBase 

♦  Model  parameters 

♦  Benchmark  results 


Other  utilities 

♦  Basic  math  calculation 

♦  File  input/output  routines 

♦  Data  format  conversion 

♦  Interface  routines 


MODELING  ROUTINES 


♦  Inverse  simulation  routines 

♦  Forward  simulation  routines 


BENCH  MODELS 


♦  Data  Analysis/interpretation 

♦  High  level  graphics 


USER  MODELS 


♦  Dynamics  bench  models 

♦  Muscle  modeling  benches 

♦  Control  benches 


APPLICATION 
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CHAPTER  2 

MULTIBODY  KINEMATICS 


Kinematics  is  the  science  of  motion  without  considering  its  relationship  with  the  force 
applied.  In  human  movement,  it  is  the  study  of  the  positions,  angles,  velocities,  and  accelerations  of 
body  segments  and  joints  during  motion.  Despite  of  the  fact  that  many  textbooks  have  been 
published  on  human  movement  (Vaughan,  Davis  et  al.  1990;  Winter  1990;  Ozkaya  and  Nordin 
1991),  there  is  no  complete  and  rigorous  formulations  for  the  kinematics  of  human  movement.  In 
addition,  controversy  stiU  exists  among  biomechanical  community  on  the  use  of  different  reference 
frames,  segment  orientation  conventions  and  the  definition  of  joint  angles. 

This  chapter  provides  the  background  information  and  detailed  kinematics  formulations. 
The  formulations  suffice  the  needs  for  most  human  movement  analysis.  The  concept  of  tree- 
structured  systems  and  various  types  of  coordinate  frames  are  described.  Formulations  for  rotational 
kinematics  are  given.  Knowledge  of  linear  algebra  is  required  to  understand  the  formulations.  More 
information  on  kinematics  are  can  be  found  in  (Meirovitch  1970;  Paul  1981;  Kane  and  Levinson 
1985). 

2. 1  TREE  STRUCTURED  MULTIBODY  SYSTEM 

A  multibody  system  may  or  may  not  contain  closed  loops.  When  there  is  no  closed  loop  in 
the  system,  it  is  called  a  tree-structured  system.  Systems  with  closed  loop  usually  have  to  be 
converted  into  tree  structured  system  by  substitute  parts  of  closed  loops  with  kinematic  constraints. 
A  tree  structured  multibody  system  has  a  single  rigid  body  called  the  root.  The  root  is  attached  to  the 
inertial  frame  through  a  joint.  Each  rigid  body  in  the  multibody  structure  has  a  unique,  non¬ 
overlapping  path  from  the  root  to  itself. 

A  graph  can  be  created  from  a  tree-structured  multibody  system  with  the  following 
conventions.  Each  rigid  body  is  assigned  as  a  node  and  each  joint  connecting  the  bodies  is  assigned 
as  an  edge.  The  inertia  frame  is  assigned  node  0.  The  root  node  is  labeled  with  the  number  1,  and  the 
remaining  rigid  bodies  are  labeled  in  a  depth-first  manner.  The  joints  are  labeled  so  that  the  joint 
index  is  the  same  as  the  unique  rigid  body  outboard  to  the  joint.  One  example  of  a  tree-structured 
representation  of  a  human  of  13  bodies  and  13  joints  is  illustrated  in  the  following  diagram. 
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Figure  2-1.  Tree-Structured  Representation  of  a  13-segment,  13- Joint  Human 

Several  coordinate  systems,  as  shown  in  Figure  2-2,  are  commonly  used  in  the  analysis  of  a 
multibody  system.  Inertia  reference  system  as  shown  in  Figure  2-2(a)  is  usually  the  global  coordinate 
system,  which  means  the  coordinates  of  its  origin  are  zero  and  all  other  coordinate  systems  are 
specified  with  respect  to  this  system.  Each  body  segment  has  a  local  coordinate  system  (marked  with 
subscript  "L"  in  Figure  2-2(b),  which  usually  has  its  origin  at  the  segment  mass  center.  The 
orientation  of  the  local  reference  system  can  be  arbitrary.  The  principal  moment  of  inertia  axes, 
subscripted  with  "P"  in  Figure  2-2(b),  are  specified  with  respect  to  the  local  reference  system.  In 
order  to  calculate  joint  angles,  it  is  necessary  to  define  two  joint  coordinate  systems  for  a  joint;  one 
rigidly  attached  to  each  of  the  two  segments  that  are  connected  by  the  joint,  as  shown  in  Figure 
2-2(d).  The  orientations  of  the  joint  coordinate  systems  are  specified  by  the  rotation  from  the  local 
reference  systems  of  both  segments.  Once  the  two  joint  coordinate  systems  are  defined,  they  are 
fixed  in  the  corresponding  segments  and  unable  to  move  relative  to  the  segments. 
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Figure  2-2.  Reference  coordinate  systems 

(a)  Global  coordinate  system;  (b)  Local  coordinate  system  and  principal  moment  of  inertia  axes  of 
segment  1  at  the  segment's  center  of  mass;  (c)  Geometry  coordinate  system  of  segment  2  at  the  segment's 

geometric  center;  (d)  Joint  coordinate  systems 

The  geometry  of  the  outer  surface  is  sometimes  defined  for  a  rigid  body  segment  to  represent 
the  physical  appearance  of  the  segment.  This  is  necessary  for  the  visual  representation  of  the  body  as 
weU  as  used  to  detect  contact  between  different  rigid  bodies  and  the  environment.  There  is  no  direct 
association  of  the  segment  inertial  properties  and  the  shape  of  the  contact  geometry.  The  geometry 
coordinate  system  as  shown  in  Figure  2-2(c)  is  also  specified  with  respect  to  the  local  reference 
system. 

2.2  REPRESENTATION  OF  RIGID  BODY  ORIENTATION 

Let  bi,  bi,  bs  form  the  right-handed  local  frame  of  a  rigid  body  B  moving  in  the  inertia  firame 
A  represented  by  ai,  3.2,  a^;  r^  is  a  vector  fixed  in  B;  and  r^  is  the  same  vector  expressed  in  A.  The 
orientation  of  B  with  respected  to  A  can  be  represented  in  a  few  ways. 

2.2.1  Rotational  (Orientation)  Matrix 

Rotational  matrix  is  defined  as  the  rotated  frame  B  expressed  in  base  frame  A,  J.e., 

R  =  >‘R«=[b,b2b3]  (2.1) 

Therefore,  the  invariance  of  vector  requires  the  vectors  expressed  in  A  and  B  frames  be 
related  by  the  following  equation 

v^  =  RV 

The  orientation  matrk  is  orthonormal,  le.,  R"^  =  I,  where  I  is  a  unit  diagonal  matrix  of 
dimension  3.  This  introduces  six  constraints  and  allows  three  independent  coordinates  representing 
the  nine  elements  of  the  rotational  matrix. 
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Also,  since  matrix  multiplication  is  not  commutative,  special  attention  should  be  paid  to  the 
sequence  of  the  multiplication  of  the  rotational  matrix,  which  indicates  the  frame  inside  which  the 
rotation  is  performed. 

♦  Post-multiplying  rotational  matrices  ...)  indicates  subsequently  rotate  with 

respect  to  rotated  frame  Bi,  B2,  etc. 

♦  Pre-multiplying  rotational  matrices  (...R^^R^B/Ryi)  indicates  always  rotate  with  respect  to 
inertia  frame  A. 

2.2.2  Euler  Angles 
Definition 

Rotation  may  also  be  described  by  three  euler  (cardan)  angles.  However,  to  uniquely 
determine  the  three  angles,  the  sequence  of  the  rotation  must  also  be  specified.  ZXZ  and  ZYX 
conventions  are  commonly  used.  In  ZXZ-convention,  as  shown  in  Figure  2-3(a),  rotations  are 
sequentially  performed  around  the  Z  X  (0,  and  Z  {\ii)  axes  respectively.  While  in  ZYX- 
convention,  as  shown  in  Figure  2-3Cb),  rotations  are  sequentially  performed  around  the  Z  ((^),  Y  (6), 
and  X  axes  respectively. 


(a)  ZXZ  Convention  Convention 

Figure  2-3.  Euler  Angles 


Conversion  between  Rotational  Matrix  and  Euler  Angles 

Let  Ci=cos(<z>),  5i=sin(<z>),  C2=cos(^,  52=sin(0,  ci=cos{\ii),  and  *=sin(^«/).  Also  designate  R^  R<) 
and  R^  to  indicate  the  rotational  matrices  around  the  individual  axis  respectively.  For  ZXZ 
convention,  rotational  matrix  can  be  obtained  by  post-multiplying  R^  R^  and  R^as 
'1  0  0 1  rc2  -^2  0]  [1  0  O' 

R^=  0  c,  -5,  ;  R^=  ^2  ^2  0  ;  ^  ^3  “-^3  (2-2) 

0  5,  c,  0  0  1  0  53  Cj 
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(2.3) 


CjCj  S.^C2S-^ 

R  =  R^(R^=  5,C3+C,C25'3  -S^S^+Cf^Ci  -C\S2 

*^2^3  ^2 

For  ZYX  convention,  the  rotational  matrix  is 

10  0  Cj  0  Sj  Cj  -$2  0 

R^=  0  c,  -5,  ;  R^=  0  1  0  R^=  5^  Cj  0  (2.4) 

0  5,  C,  -52  0  C2  0  0  1 

^1^2  ~'^1^3 '*’^1'^2‘^3  ‘^l‘^3  '*’^1‘^2^3 

R  =  R,}R^,^=  5,C2  c, 03+5,5253  -0,53+5,5203  (2.5) 

*^2  ^2*^3  ^2^3 

On  the  other  hand,  euler  angles  can  be  estimated  from  orientation  matrix.  For  ZXZ 
convention,  the  algorithm  is  as  follows 


(2.6) 


where  arctan2  is  the  extension  of  regular  arctan  function  allowing  the  calculation  of  an  angle  from  -tc 
to  Tt.  The  first  argument  of  the  function  should  be  the  sin-component  of  the  angle,  while  the  second 
component  be  the  cos-component  of  the  angle. 

For  ZYX  convention,  euler  angles  can  be  calculated  from  rotational  matrix  as 
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when  R3,  >  £•  or  Rj^  >  e 

</i  =  arctan2(R32,-R33) 

^&[-7t,7t] 

< 

6  =  arctan  2(-R3, ,  S3R32  +  C3R33 ) 

9&[-7tl2,nl2] 

\j/  —  arctan  2(s3R22  +  C3R23 ,  S3RJ2  +  ^3^23) 

(j)  =  arctan  2(-R32,R33) 

< 

6  =  arctan  2(-R3, ,  S3R32  -i-  C3R33 ) 

[-;r,-;r/2]  and  6  &  \7t ! 2^7t\ 

y/  —  arctan  2(s3R22  +  C3R23 ,  S3Rj2  +  C3R23) 

when  Rj,  <  s  and  <  s 

0 

6  -  arctan  2(-R3  j ,  R33 ) 

6 

ij/  =  arctan  2(R23,R]3) 

It  should  be  noted  that  using  euler  angles  doesn’t  limit  the  range  of  the  orientation  angles, 
thus  can  handle  tumbling  motions.  However,  using  rotational  matrix  cannot  by  itself  account  for  the 
tumbling  motion.  On  the  other  hand,  Euler  angles  suffers  from  the  so-called  “gimbal  locking”,  when 
two  of  the  three  rotational  axes  coincide  and  the  derivatives  of  the  euler  angels  are  indeterminate. 

2.2.3  Euler  Parameters 


Definition 


In  order  to  avoid  the  locking  of  euler  angles,  rotation  can  be  represented  by  quaternion 
which  uses  four  coordinates  rather  than  three.  The  basic  idea  of  the  quaternion  representation  of 
rotation  is  that  any  rotation  can  be  uniquely  defined  as  a  rotation  (p)  around  a  unit  vector  (X,).  The 
euler  parameters,  defined  as  follows,  is  one  of  the  many  quaternions 


e,  =  \  sin(|) 
=  X.2  sin(|) 
^3  =  X3  sin  A 


64  =cos(|) 


(2.8) 


Notice  since  X,  is  a  unit  vector,  the  norm  of  euler  parameters  [ei,  62,  es,  e4]^  is  one.  As  can  be  seen 
from  the  definition,  the  tumbling  motion  cannot  be  accounted  for  by  using  euler  parameters 


Conversion  between  Rotational  Matrix  and  Euler  Parameters 

Rotational  matrix  can  be  calculated  from  euler  parameters  as 
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1  -  l{e^  +  )  2(61^2  -  ^3^4 )  2(e, 63  +  ^^^4 ) 

R=  2(e,e2 +^364)  l-2(g,^+e3^)  2(62^3 -^1^4)  (2.9) 

2(61^3— ^2  ^4)  2(63^3  ^1^4)  1~2(S]  +^2  ) 

Euler  parameters  can  be  determined  from  a  given  rotational  matrix  by  the  following  scheme 


(2.10) 

2.3  ANGULAR  VELOCITY  AND  ACCELERATION 

The  angular  velocity  of  a  rigid  body  B  (bi,  b2,  bs)  moving  in  the  inertia  frame  A  (ai,  a2,  as)  is 
defined  as  follows 

=  bi — b2  bs  +b2 — ba  bi  +b3 — bi  •b2  (2. 1 1) 

dt  dt  dt 

Angular  acceleration  is  defined  as  the  time  derivative  of  the  angular  acceleration 

(2.12) 

dt  ^ 

Relationship  between  angular  velocity  of  a  rigid  body,  and  the  time  derivatives  of  euler 

angles,  <|)  =  [<{)  9  vj/]  ,  can  be  determined  from  the  invariance  of  the  angular  velocity  vector  with 
respect  to  the  reference  frame.  For  ZXZ  convention,  it  is 
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(2.13) 


For  ZYX  convention  the  relationship  is 

■-S2  0  llff 

^(0^=  S(j>=  CjSj  C3  0  -i  0  »  (2.14) 

C2C3  -^3  0  vj; 

The  inverse  relationship  is 

(2.15) 

It  must  be  noticed  that  the  T  matrix  (the  inverse  of  S  matrix)  may  be  singular  in  certain  cases. 
Physically,  the  singular  case  corresponds  to  the  “locking”  of  euler  angles,  Avhen  angular  velocity  can 
not  be  determined  from  euler  angles. 


Similarly,  the  time  derivatives  of  euler  parameters  can  be  calculated  from  segment  angular 
velocity  from  the  following  equation 


(2.16) 


2.4  TRACKING  TUMBLING  MOTION 


The  formulations  given  above  are  capable  of  handling  translation  between  the  different 
representations  of  rotation,  i.e.,  rotational  matrix,  euler  angles  or  euler  parameters.  However,  they 
are  not  exactly  the  same.  Euler  angles  allow  for  the  tumbling  motion  (where  rotation  is  not  limited 
to  one  cycle),  but  they  suffer  from  the  problem  of  “locking”.  Euler  parameters  have  no  problem  with 
locking,  but  cannot  handle  tumbling  motion  naturally. 

In  calculation  angular  velocity,  the  time  derivatives  of  either  euler  angles  or  Euler  parameters 
must  be  obtained.  Therefore  their  continuity  in  time  must  be  considered.  To  simplify  the  problem, 
rotation  is  represented  by  a  rotational  matrix.  Euler  angles  are  estimated  from  rotational  matrix, 
which  further  will  be  differentiated  and  used  to  calculate  angular  velocity.  As  can  be  seen  from 
equations  (2.6)  and  (2.7),  The  range  ^  and  ^are  limited  to  2;r  (one  circle),  whUe  the  range  of  ^is 
limited  to  n  (half  a  circle).  To  guarantee  the  continuity  of  Euler  angles,  two  extensions  to  the 
formulations  must  be  made.  First,  the  range  of  0has  to  be  extended  to  the  full  circle.  Equations  (2.6) 
and  (2.7)  give  two  set  of  formulations  to  calculate  Euler  angles  corresponding  to  cases  when  6  is 
inside  different  half  circles.  Using  formulations  for  one  half  cycle,  sudden  jumps  of  ;r  or  -;r  in  the 
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values  of  cj)  and  {/indicates  that  ^has  moved  to  the  other  half  circle.  Second,  jumps  of  27tm  Euler 
angles  have  to  be  removed.  This  requires  the  time  history  of  the  rotation  be  recorded. 

2.5  LINEAR  POSITION,  VELOCITY  AND  ACCELERATION 

Six  independent  parameters  are  required  to  describe  a  local  frame  B(bi,  hz,  bs)  with  respect  to 
the  inertial  frame  A(ai,  aa,  as),  i.e.,  three  independent  parameters  desaibing  the  orientation  of  B,  as 
described  in  the  previous  sections,  and  three  coordinates  of  the  origin  of  B  expressed  in  A.  Let  O^ 
the  origin  of  loccil  frame  B;  p^  the  position  vector  of  a  fixed  point  on  B;  and  the  same  position 
vector  expressed  in  A.  The  following  relationship  holds 

j/  =  ^R^pJ»+0^  (2.17) 

The  velocity,  v^,  and  the  acceleration,  a^,  of  the  point  expressed  in  the  inertia  frame  are 
defined  as 

v^  =  —  (2.18) 

dt 

(2.19) 

dt 

The  following  equations  relate  the  velocities  and  the  accelerations  expressed  in  A  and  B 

frames 

v^  =  v®  +  ^a)^xr®  (2.20) 

a'^  =  a^  +  '^co^  X  xr^  +  xr^  (2.21) 

where  '‘m®  and  are  the  angular  velocities  and  accelerations  of  the  local  B  frame. 

2.6  CONSTRAINTS 

The  motion  of  a  multibody  system  may  involve  two  types  of  constraints.  Configuration 
constraints  restrict  the  positions  of  each  rigid  body.  The  equations  expressing  the  restriction  are 
called  holonomic  constraint  equations.  Motion  constraints  impose  restrictions  on  the  velocities  of 
the  rigid  bodies.  When  there  is  no  motion  constraint  involved  in  a  multibody  system,  it  is  called  a 
holonomic  system;  otherwise,  it  is  a  nonholomonic  system.  The  first  and  second  order  partial 
derivatives  of  rotational  matrix,  R,  with  respect  to  euler  angles  or  euler  parameters  must  be  provided 
for  the  formulations  of  a  holomonic  system.  The  first  order  partial  derivatives  of  T  matrix  are  also 
needed  for  the  nonholonomic  systems.  There  relationships  are  given  in  the  following  tables. 

Table  2-1  1“  order  derivative  of  R  with  respect  to  ZXZ  euler  angles 


d/d^ 

d/de 

d/dY 

2^11 

^S\S%  -C\C2Ci 

SiSzSs 

-CiSi-SiCiCz 

1^21 

C\Ci  ~S\C2 

-C\S2Si 

-SiSz"f-C\C2Cz 

0 

CzSi 

SiCz 
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-^12 

S  Si  -C1C2C3 

SIS2CI 

'C\  Cl  'i~S\  C2S1 

Rii 

-CiSi  -S1C2C1 

-CIS2CI 

''SiCi-C\C2Si 

i?23 

0 

CiCi 

-S2SI 

CiSi 

SlC2 

Rz2 

S).S2 

-C1C2 

0 

i?33 

0 

0 
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Table  2-2  2"'‘  order  derivative  of  R  vv^ith  respect  to  ZXZ  euler  angles 


^/df 

^/de 

^/dy/ 

^/d(pd6 

/d^dy/ 

^/dddyr 

Rn 

-C\Cz  "^51^253 

S\ 

-C1C3  -f~SiC2S^ 

C1S2SZ 

5153-C1C2C3 

S\S2Ci 

J^21 

-SiCz‘CiC2S3 

•C\C2Si 

-Si  C3  -C\  C2S^ 

S1S2S3 

-C\S^-SiC2Cz 

'•C\S2C^ 

0 

-SiSi 

0 

0 

C2C3 

^12 

C\Sz  ~^S\C2Ci 

CiSi+SiCiCi 

^1 

S1C3  ■/■C1C253 

‘SiSzSs 

R22 

SiS^-CiCiC^ 

SiS^-CiC2Cz 

^CiCz  *^51^253 

CiS2S^ 

R23 

0 

-SiSi 

0 

0 

-C2S3 

Rzi 

-S1S2 

-SiSz 

0 

C1C2 

0 

0 

R32 

C1S2 

0 

SiCi 

0 

0 

i?33 

0 

-Cl 

0 

0 

0 

0 

Table  2-3  1”  order  derivative  of  R  ^vith  respect  to  ZYZ  euler  angles 


d/d(j> 

d/dd 

d/dyf 

Rii 

-S1C2 

-Ci&i 

0 

Rii 

C1C2 

-S1S2 

0 

Rsi 

0 

<2 

0 

Ri2 

CiC2Si 

S1S3  +C1S2C3 

R22 

SiC2Sz 

“(7i53  -^SiS2Cz 

R23 

0 

SiSz 

CiCz 

Rzi 

C1SS-S1S2CZ 

C1C2C3 

S1C3-C1S2SZ 

R32 

SiSz  '^C\S2Cz 

51C2C3 

-C1CS-S1S2SZ 

i?33 

0 

-SiCz 

CiSz 
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Table  2-4  2”“’  order  derivative  of  R  with  respect  to  ZYZ  euler  angles 


^/d(pd0 

^/d(l>dyf 

^/dOdy/ 

Ru 

-C1C2 

C1C2 

0 

S1S2 

0 

0 

Rji 

-S1C2 

-SIC2 

0 

-ClSz 

0 

0 

Rsi 

0 

Si 

0 

0 

0 

0 

Ri2 

S\C‘^-C\S2Si 

-Cl  5253 

C\S3-S\S2C3 

C1C2C3 

Ri2 

~C\C2~S\S2Si 

-5i^53 

S\S3  +C\S2C3 

51C2C3 

-^23 

0 

-C2^ 

-ClSi 

0 

0 

-S1C3 

Rz\ 

‘‘S\S2‘‘C\S2Cz 

-CxSiCi 

-5i52-Ci52C3 

-S1C2C3 

C1C3  '^S\S2S3 

-CxCiSi 

Ri2 

C\Si'S\S2Ci 

~S\S2Cz 

C1C2C3 

S\C3-C\S2S3 

‘SIC2S3 

-^33 

0 

C2C3 

-C2C3 

0 

0 

S2S3 

Table  2-5  Derivatives  of  T  with  respect  to  ZXZ  Euler  Angles 


S2>£ 

S2<£ 

d/d(l> 

d/dO 

a^dyf 

Locking 

Txi 

0 

-C2S3/S2 

C3/S2 

Tzx 

0 

0 

-53 

T„ 

0 

Si/S2^ 

-C2C/52^ 

Tx2 

0 

-C2C3/S2 

-S3/S2 

T22 

0 

0 

-C3 

T2X 

0 

-Cj/52^ 

C2S3/S2 

Tai 

0 

0 

0 

T32 

0 

0 

0 

T33 

0 

0 

0 
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Table  2-6  Derivatives  of  T  with  respect  to  ZYX  Euler  Angles 


S2>e 

S2<Z 

d/d(p 

d/de 

d/dyf 

B1 

0 

0 

0 

0 

0 

0 

nm 

0 

0 

0 

mm 

0 

S2Si/C2 

C3/C2 

Locking 

mm 

0 

0 

-S3 

0 

S2/C2 

S2C3/C2 

mm 

0 

S2C2/C2 

S2/C2 

0 

-C3  -C3 

0 

C2/C2 

-S2S2/C2 

Table  2-7  1”  order  derivative  of  R  with  respect  to  euler  parameters 


d/de\ 

d/dcz 

d/de3 

d/de^ 

Rn 

0 

-4?2 

-4?3 

0 

Rii 

2e2 

2c\ 

2e^ 

2ei 

i?31 

26i 

-2g, 

2e, 

-262 

R\2 

2e2 

2ei 

-2c\ 

-2e2 

-4?i 

-4?3 

0 

263 

262 

2C\ 

Rn 

2c, 

2€\ 

2€2 

R32 

-20, 

263 

2e2 

-2e, 

R33 

-4e, 

-4e2 

0 

0 
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Table  2-8  2”''  order  derivative  of  R  with  respect  to  euler  parameters 


B/dei 

B/d^ 

B/Bd 

^/de^dei 

B/de^Be^ 

B/Bo^Bg^ 

S'/de^G^ 

B/de^G^ 

Rii 

0 

-4 

-4 

0 

0 

0 

0 

0 

0 

0 

Rzi 

0 

0 

0 

0 

2 

0 

0 

0 

0 

2 

Rn 

0 

0 

0 

0 

0 

2 

0 

0 

-2 

0 

Rii 

0 

0 

0 

0 

2 

0 

0 

0 

0 

-2 

Rii 

-4 

0 

-4 

0 

0 

0 

0 

0 

0 

0 

Riz 

0 

0 

0 

0 

0 

2 

2 

0 

0 

Rz\ 

0 

0 

0 

0 

! 

2 

0 

0 

2 

0 

R^i 

0 

0 

0 

0 

0 

-2 

2 

0 

0 

^33 

-4 

0 

0 

0 

0 

0 

0 

0 

0 

Table  2-9  Derivative  of  T  with  respect  to  euler  parameters 


d/dci 

d/dCz 

d/dez 

d/de^ 

mm 

0 

0 

0 

0.5 

mam 

0 

0 

0.5 

0 

mm 

0 

-0.5 

0 

0 

-0.5 

0 

0 

0 

Ti2 

0 

0 

-0.5 

0 

T22 

0 

0 

0 

0.5 

Tz2 

0.5 

0 

0 

0 

T^i 

0 

-0.5 

0 

0 

Tu 

0 

0.5 

0 

0 

-0.5 

0 

0 

0 

T33 

0 

0 

0 

0.5 

T43 

0 

0 

-0.5 

0 
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CHAPTER  3 

FORWARD  DYNAMICS  FORMULATION 


In  human  motion  simulation,  human  body  is  assumed  to  be  a  multibody  constrained  system 
of  rigid  segments  connected  at  certain  joints.  There  are  a  few  ways  of  formulating  the  system.  The 
traditional  approach  has  been  to  select  a  set  of  minimal  coordinates  and  formulate  the  equations  of 
motion  as  a  second-order  system  of  ordinary  differential  equations  (ODE),  The  ODEs  are  then 
integrated  by  standard  ODE  solvers  (Meirovitch  1970;  Kane  and  Levinson  1985).  However,  recently 
effort  has  been  made  to  formulate  the  system  in  descriptor  form,  which  use  non-minimal  sets  of 
coordinates  (Lubich,  Engstler  et  al.  1995).  Descriptor  formulation  leads  to  a  set  of  differential 
algebraic  equations  (DAE)  including  equations  of  motion  (second  order  ODEs)  and  constraint 
equations  (differential  or  algebraic  equations). 

This  chapter  develops  a  numerical  method  to  solve  the  resultant  DAEs  from  descriptor 
formulation  by  converting  the  algebraic  equations  into  differential  equations.  The  relaxation  of 
position  and  velocity  constraints  is  compensated  in  the  method  by  projecting  the  integration  results 
back  to  the  original  constraint  manifolds.  This  method  can  solve  both  open-looped  and  close-looped 
systems. 

3.1  GOVERNING  EQUATION 

The  equations  of  motion  of  a  constrained  articulate  system  can  be  expressed  in  terms  of 
position  vector  p(t) ,  velocity  vector  \(t)  and  Lagrange  multipliers  A,(t)  as 

(a)  p  =  T(t,p)v 

(b)  Mv  =  f  (t,  p,  V,  >,)  -  T(p)^  G(/,  p)^  X  (4.1) 

(c)  g(^P)  =  0 

with  prescribed  initial  conditions 

P(^o)  =  Po.v(to)  =  Vo  (4.2) 

In  equation  (4.1),  (a)  relates  the  time  derivative  of  position  to  the  velocity;  (b)  indicates  the 
conservation  of  momentum;  and  (c)  represents  the  position  constraints.  Gis  the  derivative  of  the 
constraints  with  respect  to  position  variables,  J.e., 

G((,p)  =  ^ 

dp 
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Let  and  denote  the  number  of  position  variables,  the  number  of  velocity  variables 

and  the  number  of  Lagrange  multipliers  respectively.  The  interpretation  and  dimension  of  the 
variables  in  equation  (4.1)  is  summarized  in  Table  3-1. 


Table  3-1.  Definition  and  dimension  of  variables 


Variable 

Definition 

Dimension 

P 

Position 

R”' 

V 

Velocity 

R"' 

X 

Lagrange  multipliers 

g 

Constraint 

R”" 

f 

External  force 

R"' 

M 

Mass  matrix 

T 

Velocity  description  matrix 

G 

5g/5p 

The  algebraic  constraint  equation  in  (4.1)  has  to  be  converted  into  ordinary  differential 
equations.  Differentiation  of  the  (c)  of  equation  (4.1)  with  respect  to  time  leads  to  the  following 
equation 


dt 


=  gi.t+gijPj^^ 


(4.3) 


where  index  notion  is  used  and  subscript  ,t  =  dl dt  indicates  partial  differentiation  of  time  and 
subscript  ,j  =  didpj  indicates  partial  differentiation  of  the  y""  component  of  the  position  vector. 
Equation  (4.3)  replaces  the  position  constraints  by  the  equivalent  velocity  constraints. 


Differentiating  (4.3)  again  with  respect  to  time  leads  to 


=  giM  +  gi.,Pj  +  [gij  )  ,  Pj  +  gij  (Pj  ) , 


=  0 


Notice  Pj  =  TyVj  and  T^j,  =  0 ,  equation  (4.3)  can  be  rewritten  as 

gijTj,v,  =  - {g. „  Ig^^Pj  +  g,_j,pjp,  +  gijTj,ip,v, }  (4.4) 

Equation  (4.4)  replaces  the  position  constraints  by  acceleration  constraints.  Therefore,  the 
ODE  formulation  of  the  constrained  articulated  system  takes  the  following  form 
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3.2  NUMERICAL  FORMULATION 


M  is  a  by  square  matrix.  When  formulated  in  natural  coordinate  system,  M  is  also 
block  diagonal  and  positive  definite.  Therefore,  (b)  of  equation  (4.5)  can  be  rewritten  as 

v  =  M-'{f-T^G^X} 

By  substituting  it  into  (c),  the  explicit  form  of  "k  is  obtained  as 
k  =  (GTM-‘T''G^  )■’  {GTM-'f  +  g'  +  g"  +  g‘ } 

Notice  that  GTM'^T’^G^  is  always  a  positive  definite  square  matrix.  Therefore,  the 
following  standard  system  of  first-order  equations  is  obtained 


(4.6) 


Equation  (4.6)  can  be  solved  using  the  ODE  solvers  in  Matlab. 

3.3  PROJECTION  OF  POSITION  AND  VELOCITY  CONSTRAINTS 

Equation  (4.5)  and  its  numericcd  counterpart  equation  (4.6),  the  position  constraints  (the  (c) 
of  equation  (4.1))  are  replaced  by  the  acceleration  constraints.  During  numerical  integration, 
numerical  drifting  may  occur  and  the  position  constraints  and  velocity  constraints  (4.3)  may  not  be 
accurately  satisfied.  The  projection  method  (Lubich,  Engstler  et  al.  1995)  is  used  to  project  the 
solution  to  the  position  and  velocity  constraint  manifolds.  The  projection  equations  are  given  as 
follows. 

3.3.1  Projection  of  position  constraint 

After  a  successful  integration  step,  an  approximate  position  solution  p®  is  obtained.  It  is 
projected  to  the  position  constraint  g(^p)  =  0by  performing  the  following  Newton-Ralphson 
iterations 

(4.7) 
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until  ||t(p'’  )Ay'‘  |  <  tolerance . 

3.3.2  Projection  of  velocity  constraint 

Velocity  constraint  (4.3)  can  also  be  written  in  matrix  form  as  g'  +  Gv  =  0 .  After  a 

successful  integration  step,  an  approximate  velocity  solution  v°  is  obtained.  The  velocity  constraint 
is  implemented  by  the  following  projection 
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CHAPTER  4 

INVERSE  DYNAMICS  ANALYSIS 


When  human  body  is  assumed  to  be  a  system  of  rigid  segments  connected  at  certain  joints, 
human  motion  is  governed  by  a  set  of  equations  of  motion  which  include  body  segment  parameters, 
driving  forces  and  moments,  and  kinematic  data.  In  the  previous  chapter,  the  equations  are  solved 
for  the  kinematics  given  the  body  segment  parameters  and  the  driving  forces.  The  equations  can  also 
be  solved  for  driving  forces  given  the  body  segment  parameters  and  the  kinematics.  This  is  called 
inverse  dynamics  approach. 

Inverse  dynamics  approach  is  commonly  used  in  biomechanical  analysis  of  human  motion 
to  estimate  the  mechanical  loads  on  Joint  and  segments.  Besides  rigid  body  assumption,  it  is  also 
assumed  that  no  translational  movement  can  occur  at  joint.  The  joint  forces  and  moments  calculated 
from  inverse  dynamics  analysis  are  the  resultant  forces  or  moments  from  active  skeletal  muscles  and 
passive  joint  structures  such  as  ligaments  and  joint  capsules.  The  calculated  forces  and  moments  can 
be  further  processed  to  yield  muscle  forces  as  well  as  some  energetics  quantities  such  as  internal 
work,  external  work,  and  joint  power  flow.  These  kinetics,  energetics  and  muscle  quantities  can  be 
used  to  analyze  both  the  healthy  and  pathological  gait,  assist  in  the  diagnosis  of  the  underlying 
pathologies  of  abnormal  gait  pattern,  and  help  the  evaluation  and  design  of  soldier  equipment.  This 
chapter  describes  some  common  issues  of  inverse  dynamics  analysis. 

4. 1  MEASURING  KINEMATIC  DATA  AND  GROUND  REACTION  FORCE 

The  kinematic  data  required  for  inverse  dynamic  analysis  include  the  positions,  velocities 
and  accelerations  of  segment  centers  of  gravity  (CG)  and  rotational  angles,  angular  velocities  and 
angular  accelerations  of  all  segments.  The  displacements  of  body  segments  can  be  measured  from  a 
motion  capture  system.  The  displacements  are  then  differentiated  to  derive  the  velocities  and 
accelerations  of  the  segments.  Goniometers  and  accelerometers  can  be  used  to  complement  the 
estimation  from  the  motion  capture  system  (Ladin  and  Wu  1991;  van  den  Bogert,  Read  et  al.  1996). 

Ground  reaction  forces  can  be  measured  by  force  plates,  foot  contact  pressure  measurement 
devices  (Han,  Paik  et  al.  1999)  or  alternately  estimated  from  the  kinematic  data  (Bobbert, 
Schamhardt  et  al.  1991).  The  kinematic  data  and  the  ground  reaction  data  must  be  synchronized. 

4.2  KINEMATICS  RECONSTRUCTION 

A  video  based  motion  capture  system  records  the  positions  of  clusters  of  markers  attached  on 
the  surfaces  of  body  segments  in  the  global  (laboratory)  reference  system.  The  recorded  coordinates 
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are  used  to  reconstruct  the  positions  of  anatomical  landmarks  Ooint  positions)  and  the  orientation  of 
the  segments  and  thus  the  positions  of  segment  centers  of  gravity. 

At  least  three  markers  are  needed  for  each  segment.  Given  the  position  of  three  non-collinear 
markers  on  a  segment  and  assume  that  there  is  no  relative  movement  between  the  markers,  the  joint 
position  can  easily  be  reconstructed.  The  orientation  of  the  segment  can  be  calculated  by  geometric 
rules.  Both  the  distance  between  the  markers  associated  with  each  marker  and  the  offset  of  each 
marker  form  the  line  connecting  any  other  two  should  be  sufficiently  large  to  prevent  excessive  error 
propagation  during  reconstruction.  However  since  the  markers  are  placed  on  the  surfaces  of 
segments,  relative  movement  of  markers  (skin  marker  artifacts)  may  exist.  It  has  been  shown  that  the 
errors  in  kinematic  reconstruction  due  the  skin  marker  artifacts  can  be  overwhelming  (Cappozzo, 
Cappello  et  al.  1997). 

In  order  to  remove  the  effects  of  the  skin  movement  artifacts,  a  few  techniques  have  been 
developed  by  putting  more  markers  one  each  segment.  The  redundancy  of  the  marker  data  allows 
for  an  optimal  identification  by  using  either  a  singular  value  decomposition  algorithm  or  a  weighted 
least  square  algorithm  (Spoor  and  Veldpaus  1980;  Veldpaus,  Woltring  et  al.  1988;  Soderkvist  and 
Wedin  1993). 

4.3  BODY  SEGMENT  PARAMETERS 

The  body  segment  parameters  refer  to  the  inertial  parameters  of  a  body  segment,  such  as  the 
mass  of  the  segment,  its  length,  the  location  of  its  center  of  gravity  with  respect  to  the  segment  local 
coordinate  frame  and  moments  of  inertia  (the  values  and  orientation  of  the  principal  moments  of 
inertia).  The  accuracy  of  the  body  segment  parameters  affects  the  accuracy  in  the  kinetic  analysis 
(Krabbe,  Farkas  et  al.  1997). 

A  number  of  ways  have  been  used  to  estimate  the  body  segment  parameters.  The  inertia 
properties  can  be  measured  directly  from  cadavers  (Dempster  1955;  Chandler,  Clauser  et  al.  1975). 
However,  these  data  include  only  limited  sample  size,  age  and  morphological  discrepancy.  Mass  and 
density  of  individual  segments  can  also  be  determined  from  direct  measurements  on  living  subjects. 
A  number  of  techniques  such  as  magnetic  resonance  imaging  (Martin,  Mungiole  et  al.  1989; 
Mungiole  and  Martin  1990)  and  computerized  tomography  (Brooks  and  Jacobs  1975;  Huang  and 
Wu  1976;  Zatsiorsky  and  Seluyanov  1985)  have  been  used.  These  techniques  can  provide  very 
accurate  in  vivo  estimates  of  the  segmental  inertial  properties.  However,  these  techniques  are 
difficult  to  perform  due  to  exposure  to  radiation,  high  cost  and  lengthy  procedure.  Also,  the  direct 
measurement  of  inertial  properties  can  not  satisfy  die  need  in  cissessing  all  the  body  segment 
parameters  on  individual  subjects.  Therefore,  many  linear  and  nonlinear  statistical  regression 
equations  are  developed  based  on  the  measured  data  (Hinrichs  1985;  Yeadon  and  Morlock  1989; 
Vaughan,  Davis  et  al.  1990;  Zatsiorsky,  Seluyanov  et  al.  1990).  In  addition,  body  segment 
parameters  can  be  calculated  mathematically  using  the  assumption  that  the  segment  can  be 
represented  by  simple  geometry  with  uniform  density  (Jensen  1978;  Hatze  1980;  Jensen  1989).  High 
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accuracy  has  been  reported  and  an  automated  version  is  developed  based  on  video  image  capturing 
and  computer  image  processing  (Sarfaty  and  Ladin  1993). 

4.4  JOINT  DYNAMICS  FORMULATION 


Figure  4-1  Free  diagram  of  a  segment 

Figure  4-1  shows  the  free  diagram  of  a  segment  with  global  frame,  segment  local  frame, 
external  forces  and  joint  reaction  forces  displayed.  Coordinate  system  0(X,Y,Z)  is  the  global 
reference  frame.  Local  coordinate  system  is  defined  at  the  segment  center  of  gravity  o,  as  Os(x,y,z). 
Therefore,  orientation  of  the  segment  can  be  represented  by  the  rotational  matrix  R'  =  =  [x  y  z]. 

Points  Pand  Z?/are  the  proximal  and  distal  joints  of  the  segment.  Subscript  /indicates  more  than  one 
distal  joints  may  be  present  on  one  segment  (say,  two  hip  joints  as  distal  joints  of  trunk  segment). 
However,  only  one  proximal  joint  is  allowed. 

At  the  proximal  joint,  reaction  force  vector  F^  and  reaction  torque  vector  T^  have  positive 
poles,  which  means  the  positive  direction  of  the  force  or  torque  components  is  the  same  as  the 
positive  direction  of  the  corresponding  local  coordinate.  While  at  the  distal  joint(s),  reaction  force 
vector  Fj/and  reaction  torque  vector  T^,  have  negative  poles.  Vectors  Tp  and  ta  represent  the  relative 
positions  of  the  proximal  and  distal  joints  respectively.  Force  vector  F,  is  the  additional  external 
force  acting  on  the  segment.  Several  or  none  F^  may  be  involved  depending  on  the  situation. 
Similarly,  vector  T,  stands  for  the  external  torque  exerted  on  the  segment  by  the  environment.  One 
example  of  tiie  external  force  is  the  ground  reaction  forces  acting  on  ankles.  Vector  r^  shows  the 
relative  position  of  the  application  point  of  external  force  to  o,. 

Let  m*  indicate  the  mass  of  the  segment;  L  the  inertia  matrix  of  the  segment  as  described  in 
the  previous  section;  and  g  the  gravity  vector.  Let  a*,  Oj  and  a*  be  the  linear  acceleration  vector. 
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angular  velocity  vector  and  angular  acceleration  vector  of  the  segment,  represented  in  the  global 
frame,  respectively.  The  equation  of  motion  can  be  obtained  by  conservation  of  momentum  as 
follows 

F„  +  -2  Frf/  =  m^a,  (4. 1) 

T,  +  Tp  +  Tp  xT,.  -2  TdiX  Tdi  =  Is-tts  +c)sx(Is-o)s)  (4.2) 

Since  kinematic  data  a*,  tOj  and  as  are  obtained  from  kinematic  reconstruction  and  the 
differentiation  of  the  position  data,  equations  (4.1)  and  (4.2)  for  each  segment  can  be  solved  by 
starting  with  the  most  distal  segment  and  proceeding  to  the  most  proximal  one. 

The  resultant  joint  reaction  forces  and  torques  calculated  from  equations  (4.1)  and  (4.2)  are 
defined  in  global  reference  coordinate  system,  which  are  usually  not  suitable  for  interpreting  relative 
to  the  subject’s  joints.  Therefore,  the  joint  forces  and  torques  are  usually  transformed  to  a  body  fixed 
anatomical  reference  frame.  The  anatomical  coordinate  system  first  proposed  by  [Grood,  1983  #1]  is 
widely  used  [Apkarian,  1989  #5;  Chen,  1988  #4;  Siegler,  1988  #3;  Vaughan,  1990  #6].  Also,  it 
should  be  noticed  that  the  calculated  joint  loads  are  the  resultant  loads  from  active  muscles  as  well 
as  the  passive  structures. 
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CHAPTER  5 

EXAMPLES  AND  APPLICATIONS 


This  chapter  gives  sample  applications  and  example  models  developed  from  AHBM 
V  1.  Among  inverse  dynamic  applications,  GaitSD  VI,  a  three-dimensional  gait  analysis 
program  is  developed  for  USARIEM;  a  new  method  is  developed  for  a  DOT  project  to 
understand  the  airbag  external  load  behavior  and  bag-occupant  interaction;  and  a  three- 
dimensional  lower  extremity  benchmark  model  is  developed.  Among  forward  dynamic 
examples,  a  simple  spring  mass  model  and  a  double  pendulum  model  are  developed  to  test 
the  kinematics  and  integration  algorithms  of  the  toolbox.  A  three-dimensional  whole  body 
human  model  and  a  two-dimensional  human  head-neck  model  are  developed  for  impact 
simulations. 

5.1  INVERSE  DYNAMICS  EXAMPLES  AND  APPLICATIONS 

5.1.1  GaitSD  VI  Developed  for  USAREIM 

The  kinematics,  inverse  d5mamics  and  graphical  algorithms  of  the  AHBM  toolbox 
are  used  to  develop  a  three-dimensional  gait  analysis  program  for  the  U.S.  Army  Research 
Institute  of  Environmental  Medicine  (USARIEM).  The  version-1  of  the  program,  GaitSD  VI 
(Shen  2000),  supports  the  current  test  setup  of  USARIEM  motion  analysis  system.  Input 
and  output  file  formats  currently  used  by  USARIEM  are  fully  supported. 

GaitSD  VI  allows  running,  visualizing  and  analyzing  multiple  gait  tests  in  a  simple 
graphical  layout.  It  supports  running  multiple  tests  in  one  project  and  provides  easy  control 
of  analysis  parameters.  GaitSD  VI  also  comes  with  the  StickViewer  and  XYviewer 
developed  in  the  AHBM  toolbox.  These  viewers  allow  the  visualization  and  animation  of 
test  results. 

If  requested  by  USARIEM,  the  program  will  be  upgraded  to  include  the  further 
development  of  AHBM  toolbox  in  the  following  phases. 


SS 


(a)  Graphical  Layout  to  Edit  Multiple  Gait  Tests 


GAIT  PROGRAM;  |D;\Gor^Piogiam\Vetsionl \Exomple\T cslPiojeclVsamplu.proi] 


pErojed  led  J$<V5ewer 

Project! 


(b)  StickerViewer  and  XYViewer 


Figure  5-1.  GaitSD  VI  Developed  for  USARIEM 
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5.1.2  Inverse  Dynamics  Model  for  Predicting  Airbag  Forces 

An  inverse  model  as  shown  in  Figure  5-2,  is  developed  for  calculating  external  air 
bag  loads  on  the  head  and  neck  of  a  small  female  test  dummy  using  recorded  dummy 
response  data  (Chan,  Shen  et  al.  2001).  This  work  is  sponsored  by  Department  of 
Transportation  to  understand  and  minimize  the  possible  hazard  of  air  bags  based  on 
physical  principles. 


Figure  5-2.  Schematics  of  the  Inverse  Model  to  Predict  Airbag  Force 

Calculations  are  performed  for  static  out-of-position  tests  as  well  as  vehicle  crash 
tests.  Some  sample  calculations  are  given  in  Figure  5-3.  The  calculated  external  loads 
provide  a  phenomenological  explanation  of  the  differences  in  dummy  responses  between 
different  dummy  positions.  It  is  also  shown  that  the  impact  angle  on  the  head  affects  the 
head/neck  joint  load  significantly. 


(a).  Head  Rotation  Angle  (deg)  (b).  Airbag  Load  on  Head  (KN) 

Figure  5-3.  Calculated  Head  Rotation  and  Airbag  Load 
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5.1.3  A  three-dimensional  lower  extremity  model 

A  three-dimensional  lower  extremity  model  (Vaughan,  Davis  et  al.  1990)  is 
developed.  This  model  consists  of  seven  segments,  i.e.,  feet,  shanks,  thighs  and  a  pelvis.  Six 
joints,  ankles,  knees  and  hips,  connect  the  seven  segments.  In  the  bench  problem  described 
in  (Vaughan,  Davis  et  al.  1990),  the  three-dimensional  motion  (positions  and  orientations) 
of  the  lower  extremity  is  reconstructed  from  eighteen  markers  attached  to  the  segments. 
Two  force  plates  record  the  ground  reaction  forces  on  both  feet.  The  mass  and  inertia 
properties  of  each  segment  are  calculated  from  regression  equations. 


Figure  5-4.  A  three-dimensional  lower  extremity  model 

The  kinematics  and  ground  reaction  forces  given  in  (Vaughan,  Davis  et  al.  1990)  are 
used  to  test  the  kinematic  and  inverse  d3mamic  routines  in  the  toolbox.  Results  are 
consistent  to  those  given  in  the  book.  Some  of  the  comparison  of  result  is  given  in  Figure 
5-5. 

This  model  can  be  used  to  perform  three-dimensional  gait  analysis.  Lower  extremity 
muscles  will  be  added  to  the  model  in  the  next  phase  of  AHBM  development.  This  will  allow 
the  comparison  and  testing  of  different  muscle  models  and  muscle  functions  during  human 
locomotion. 
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(a)  Proximal/distal  force  at  the  right  ankle  joint  in  Newton.  Red  line  is  the  result  calculated 
using  toolbox  routines;  blue  line  is  from  (Vaughan,  Davis  et  al.  1990) 


(b)  Flexion/Extension  torque  at  the  right  ankle  joint  in  (Nm).  Red  line  is  the  result 
calculated  using  toolbox  routines;  blue  line  is  from  (Vaughan,  Davis  et  al.  1990) 

Figure  5-5.  Comparison  of  Forces  and  Torques  at  the  Right  ankle 
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5.2  FORWARD  DYNAMICS  ANALYSIS  EXAMPLES 


5.2.1  Simple  Spring  Mass  Model 

Figure  5-6  shows  a  simple  model  where  two  masses,  mj  and  m^,  are  connect  hy  a 
linear  spring  with  a  spring  coefficient  of  k.  A  time  dependent  force  F(t)  is  applied  on  mass 
mj.  This  problem  is  used  to  test  the  numerical  algorithm  in  the  toolbox  to  solve  equations  of 
motion.  Figure  5-7(b)-(d)  give  the  solution  under  a  cyclic  load  F(t)  =  10sin(2*TC  *t)  (Figure 
5-7(a))  and  mj  =  10  (kg);  mj  =  10  (kg);  k  =  10  (N/m). 


m  M  St  jm  Mr  m  m  m  m 


Figure  5*6.  A  mass-spring-mass  system 


(a).  External  Force 


0  0.5  1  1.5  2 


(b).  Spring  reaction  force 


Time  (second) 

(d).  Position  of  Mass  #2 


Figure  5-7.  External  force  and  solutions 
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5.2.2  Three-dimensional  Double  Pendulums 


A  three-dimensional  double  pendulum  is  shown  in  Figure  5-8.  The  lower  pendulum 
is  connected  to  the  ground  by  a  pin  joint,  where  rotation  is  allowed  around  only  one  axis  of 
the  joint.  An  Euler  joint  connects  the  upper  and  the  lower  pendulum  and  allows  for 
subsequent  rotation  around  all  three  axes  of  the  joint.  This 
problem  is  used  to  test  the  algorithms  in  the  toolbox  to 
calculate  three-dimensional  kinematics. 

In  this  example,  both  pendulums  are  of  unit 
lengths,  imit  masses  and  xuiit  moments  of  inertia.  A 
constant  unit  force  is  applied  at  the  top  of  the  upper 
pendulum  and  moves  with  the  pendulum  (expressed  in  the 
upper  pendulum  local  frame).  The  calculated  kinematics  is 
given  in  Figure  5-9(a)-(d). 

Figure  5-8.  3D  Double  Pendulum 


0  0.2  0.4  0.6  0.8  1 


(a)  Lower  (pin)  joint  angle  (deg)  (b)  Upper  pendulum  ang.  vel.  (rad/s) 


0  0.2  0.4  0.6  0.8  1 

Time  (second) 


(c)  Upper  (Euler)  joint  angles  (deg) 


Time  (second) 

(d)  Upper  pendulum  vel.  (m/s) 


Figure  5-9.  Calculated  Kinematics  of  the  Double  Pendulum 
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5.2.3  Three-dimensional  Whole  Body  Human  Model 

Figure  5-10  shows  a  three-dimensional  whole  body  model  developed  using  AHBM  VI 
routines.  The  model  has  15  body  segments,  i.e.,  head,  neck,  upper  torso  (thorax),  center 
torso  (abdomen),  lower  torso  (pelvis),  upper  arms,  lower  arms,  upper  legs,  lower  legs  and 
feet.  The  lower  arms  are  combinations  of  the  forearms  and  hands.  The  are  connected 
together  by  14  joints  representing  the  physical  joints  of  the  human  body  such  as  pelvis, 
waist,  neck,  hips,  knees,  ankles,  shoulders  and  elbows. 

Various  types  of  joints,  such  as  pin  joint,  ball  and  socket  joint  and  euler  joint,  are 
used.  The  joints  are  constrained  by  nonlinear  rotational  springs  and  dampers.  Joint  range 
of  motion  is  implemented  by  including  “soft  joint  stop”,  i.e.,  when  angles  beyond  the  joint 
range  of  motion,  a  restoring  torque  of  significant  magnitude  is  added  to  move  the  joint  back 
to  the  range  of  motion. 


Figure  5-10.  Three-dimensional  15-segment  human  model 

The  body  segment  inertia  and  geometrical  properties  and  joint  mechanical 
properties  refer  to  (Cheng,  Obergefell  et  al.  1994).  These  values  are  easily  adjustable  due  to 
the  flexibility  of  toolbox  routines.  The  example  used  here  has  a  weight  of  about  75 
kilograms  and  a  height  of  about  1.8  meters.  Based  on  these  two  parameters,  it  is  classified 
in  the  95*'’  percentile  of  the  adult  U.S.  males 

An  impact  force  as  shown  in  Figure  5-ll(a)  is  applied  at  the  chest.  The  model  is  used 
to  calculate  the  response  of  the  whole  body  to  the  impact.  The  chest  acceleration,  head 
acceleration  and  head  angular  accelerations  are  given  in  Figure  5-ll(b)-(d). 
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(a)  Impact  Force  on  Chest  (KN) 


(b)  Chest  Acceleration  (g) 


Time  (millisecond)  Time  (millisecond) 

(c)  Head  Angular  Acceleration  (rad/s^)  (d)  Head  Acceleration  (g) 

Figure  5-11.  Impact  Force  and  Calculated  Accelerations 

This  model  shows  the  capability  of  toolbox  routines  in  dealing  with  impact  forward 
simulation  problems  of  large  degrees  of  freedom  and  complex  joint  constraints. 
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5.2.4  Two-dimensional  Human  Head-Neck  Model 

A  human  head-neck  model  in  is  developed  using  AHBM  toolbox  routines.  This  model 
is  a  simplified  version  of  the  global  human  head-neck  model  described  by  (De  Jager  1996). 
As  shown  in  Figure  5-12,  nine  rigid  bodies  represent  the  head  (CO),  the  seven  cervical 
vertebrae  (C1-C7)  and  the  first  thoracic  vertebra  (Tl).  The  inertia  and  geometrical 
properties  of  each  segment  were  measured  from  test  and  given  in  (De  Jager  1996).  The 
bodies  are  connected  by  complex  nonlinear  viscoelastic  elements. 


Figure  5-12.  Two-dimensional  Human  Head  Neck  Model 

The  model  is  used  to  simulate  the  head-neck  response  during  a  front  impact  test. 
The  sled  acceleration  history  is  given  in  Figure  5-13(a).  The  simulation  results  are 
consistent  with  those  given  in  (De  Jager  1996).  Part  of  the  simulation  results,  such  as  total 
head  rotation  angle  and  joint  angles  between  some  neighboring  cervical  vertebra  are  given 
in  Figure  5-13(b)-(f). 

Since  the  neck  joints  are  connected  by  nonlinear  springs  which  may  become  very 
stiff  at  certain  alignment  and  orientation.  This  model  demonstrates  the  toolbox’s  ability  of 
solving  stiff  problems. 

The  simulation  results  also  show  that  the  effects  of  active  neck  muscles  must  be 
included  for  accurate  prediction  of  head-neck  response,  even  under  very  violent  situations 
when  the  muscle  effects  are  usually  neglected. 
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(a)  Sled  Acceleration  (on  Tl)  (m/s^)  (b)  Head  Rotation  Angle  (deg) 


(c)  Joint  angle  between  C7-T1  (deg) 


(d)  Joint  angle  between  C6-C7  (deg) 


(e)  Joint  angle  between  C5-C6  (deg)  (f)  Joint  angle  between  C4-C5  (deg) 


Figure  5-13.  Sled  Acceleration,  Head  and  Joint  Angles 
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